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Abstract. We discuss work toward developing a 2.5D non-LTE radiative 
transfer code. Our code uses the short characteristic method with modifications 
to handle realistic 3D wind velocities. We also designed this code for parallel 
computing facilities by representing the characteristics with an impact parameter 
vector p. This makes the interpolation in the radiation angles unnecessary and 
allows for an independent calculation for characterisitcs represented by different 
p vectors. The effects of the velocity field are allowed for by increasing, as 
needed, the number of grid points along a short characteristic. This allows us 
to accurately map the variation of the opacities and emissivities as a function 
of frequency and spatial coordinates. In the future we plan to use this tran sfer 
code with a 2D non-LTE stellar atmosphere program l)Georgiev et al.ll2004(l to 
sclf-consistently solve for level populations, the radiation field and temperature 
structure for stars with winds and without spherical symmetry. 



1. Introduction 

Modeling the circumstellar envelopes of O and B stars is a complex nonlinear 
problem. The non-LTE level populations, the (magneto-) hydrodynamics, and 
the radiation field are strongly coupled. Thus, an iterative procedure is needed 
to achieve a consistent solution. An essential constituent of this procedure is 
the availability of an accurate and fast radiative transfer code. 

Progress in computer technology and the availability of fast numerical meth- 
ods now allow the development of such codes for detailed study of 2D and 3D 
envelopes. There are several possible avenues to follow. The most straightfor- 
ward is to solve the general radiative transfer equation 



and calculate the radiative transition rates using the solution. In the above n 
is the direction of the radiation, I is the specific intensity, x 1S the opacity, and 
S is the source function (functional dependence is not indicated for clarity). 
Alternatively, one could also use the moment equation s, derived from Eq. and 
solve directly for the moments |Hilher &; Miller] Il998h which set the transition 
rates; or use Monte-Carlo simulation to solve the radiative transfer equation 
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and calculate estimators of the transition rates (|LucvHl999T ). We decided to use 
the first approach because of its simplicity and since it provides a reasonable 
compromise between numerical efficiency and flexibility. 

A simple iteration between the radiative transfer and the rate equations is 
not a wise choice for the iterative procedure. This is the so-called "A-iteration" 
which is notorious for convergence difficulties when the optical depth is large. 
Convergence i s ensured, however, by using the Approx imate Lambda Iteration 
(ALI, see e.g., iRvbicki fc Hummer! Il99ll : lHubenvHl992T l which takes some cou- 
pling of the radiation and populations into account by using an invertible Ap- 
proximate Lambda Operator (ALO). In our 2D code we use the local contribu- 
tion at every spatial point to construct the ALO because it is easy to calculate 
and has acceptable convergence characteristics. The actual implementation of 
the ALI procedure int o our full non-LTE model atmosphere will be discussed in 
iGeorgiev etafl (|2004T l. 



2. Description of the Transfer Code 

The optical depth and the formal solution of Eq. ^at any position, s, along a 
ray are 

Tu = f Xds' and I(t u ) = I BC e~ Tv + [" S(r') e T '~ T » dr' , (2) 
Jo Jo 

respectively. The intensity can be calculated by specifying Ibc at the up-stream 
end of the ray (or characterisitic) and by evaluating two integrals. For this 
purpose, we use the "Sh ort Characteristi c " (SC ) method, first explored by 
iMihalas etaD (|l978h and iKunasz & Auerl fl988). This method requires the 
evaluation of the integrals only between the point of interest and the closest cell 
boundary and uses the calculated intensities at other grid points to interpolate 
Ibc ■ in the spherical coordinate system the directional variation of the intensity 
is normally described by the radiation coordinates 9 and </>, which are defined 
by 

cos(9) = n • r and sin(/3) ■ sin(9) ■ cos{4>) = [n x r ] • [r x z] , (3) 

respectively (see Fig. ^ for definitions). Unfortunately, 9 and eft vary along a 
characteristic so using the same 9 and <j) grid for all spatial points would require 
interpolations in these angles. To avoid this additional interpolation we describe 
a characteristic with 

p = r x n , (4) 

which we call the "impact parameter vector" (see Fig. ^). This vector describes 
all essential features of a characteristic and can be considered as an analog of 
the orbital momentum vector in two body problems. Its absolute value p= 
|p| is the traditional impact parameter while its orientation defines the "orbital 
plane" of the radiation (the plane that contains the characteristic and the origin) . 
Following the analogy one can define an "inclination" angle for this plane by 

p ■ cos(i) = p • z . (5) 
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In our code we set up a universal grid in impact parameters and in inclination 
angles and calculate the radiation coordinates by 



for each grid point (see Fig ^ for definitions). 

We evaluate Eq. [2] in the comoving frame of the point of interest which is 
the proper frame for solving the rate equations. To ensure that the spatial and 
frequency variations of the opacity and source function are mapped properly 
in the integrations, we add extra integration points to the characteristics. The 
number of the extra points (at least one) depends on the ratio of the line of 
sight velocity difference between the endpoints and a "maximum allowed velocity 
difference" which is a free parameter in the code. The opacities and source terms 
at every comoving frequency are then interpoleted onto the integration points by 
bi-linear approximations using the four closest spatial grid points. It is relatively 
easy to construct a diagonal ALO in this evaluation scheme. 

With the exception of the intensity we interpolate all quantities in first 
order. However, the accuracy of this approximation is insufficient in many cases; 
therefore, we introduced a rudimentary multi-grid approach. Before evaluating 
Eq. [2J we calculate opacity and source function on a dense grid by using a 
sophisticated 3 rd order approximation. Then, the transfer equation is solved on 
the original grid using the dense grid for opacity and source term interpolations. 

3. Test Results 

We te sted our code by reproducing spherical symmetric CMFGEN models ((Hillier Sz Millerl 
199S [), as well as the re sults of an accurate 2D long characteristic program (see 
Bus che fc HiilieU I2000T ) . The 2D models were static with Schuster-type inner 
boundary conditions and included electron scattering iterations. We were able 
to reproduce all test cases within ~5% accuracy. In Figs. |2 and |3] we demon- 
strate the capabilities of our code by showing the results for a rotating stellar 
envelope. The model was produced by using the opacities and emissivities from 
the results of a realistic CMFGEN simulation and by introducing a rotational 
velocity field. As expected the spectral lines show the rotational broadening. 
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Figure 1. The definition of our fundamental coordinate system (r, /3, e) 
and the vectors that are necessary for our notations. Unit vectors n, z, and r 
are describing the characteristic (long thin line), the positive z axis, and the 
radial direction, respectively. The impact parameter vector p and n are two 
alternative ways to define a characteristic. 
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Figure 2. The J moment of the radiation field for j3= 0° (thin line) and 
45° (thick line) as a function of wavelenght (A) at the outer boundary of the 
envelope. 
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Figure 3. The spectral region around 1444 A from Fig. El as a function 
of velocity (centered on 1444 A). The rotational velocity in the envelope is 
V(r, 0) = 250 • r*/r ■ sin{J3) where r* is the stellar radius. 



